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Bilayer graphene is a recently isolated and intriguing class of many-body systems with massive 
chiral quasiparticles. We present theoretical results for the electronic compressibility of bilayer 
graphene that are based on a four-band continuum band structure model combined with a ran- 
dom phase approximation treatment of electronic correlations. We find that the compressibility is 
strongly suppressed by electron-electron interactions at low carrier densities. Correlations do not 
lead to any qualitative new features, but are crucially important for a quantitative understanding 
of this fundamental thermodynamic property of graphene bilayers. 
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I. INTRODUCTION 

Crystalline bilayers of graphene (BLG) produced by 
mechanical exfoliation of thin graphite or by thermal 
decomposition of silicon carbide have recently attracted 
a great deal of attention because of their many unique 
electronic properties^"*. BLG quasiparticles behave at 
low energies like massive chiral fermions^ and exhibit 
a plethora of interesting properties, including broken- 
symmetry states at very vi^eak magnetic fields when the 
bilayer is suspended^ to reduce disorder, and anomalous 
exciton condensation in the quantum Hall regime^°. 

Since BLG consists of two single-layer graphene (SLG) 
systems separated by a small distance d ~ 3.35 A, one 
expects inter-layer electron-electron interactions be cru- 
cial to the physics of this system. With this motivation, 
many-body effects in BLG have already been studied by 
several authors^^"^°. Particular attention has been de- 
voted to the study of interaction effects close to charge 
neutrality^-'^"^^ (see also Ref. 26) where it has been shown 
that BLG is prone to a number of interesting instabili- 
ties, including sublattice pseudospin ferromagnetism - a 
type of orbital order which leads to spontaneous inversion 
symmetry^"'^'^^ breaking. 

Thermodynamic quantities such as the electronic com- 
pressibility K or the spin susceptibility xs are very pow- 
erful probes of exchange and correlation effects in in- 
teracting many-electron systems^^ since they are inti- 
mately linked with the equation of state. The elec- 
tronic compressibility of a conventional parabolic-band 
two-dimensional (2D) electron gas was first measured by 
Eisenstein et al?^ in 1992. For sufficiently low densities, 
and zero magnetic fields, it was found that the inverse 
thermodynamic density-of-states, which is proportional 
to 1/k, changes sign becoming negative, a fact that can 
be easily explained by properly including exchange con- 
tributions to the free-electron equation of state^^. In an 
ordinary 2D electron gas corrections to the compress- 
ibility due to correlation effects omitted in Hartree-Fock 
approximations are relatively small. The "field penetra- 



tion technique" introduced in Ref. 28 and later discussed 
in great detail in Ref. 29 actually uses a double-layer 2D 
electron system made up of two closely-spaced 2D elec- 
tron gases and can also be used to accurately measure 
the compressibility of BLG. 

In this work we present a calculation of the electronic 
compressibility of BLG based on the four-band contin- 
uum model. We include beyond-Hartree-Fock correla- 
tion contributions to the ground-state energy by using a 
random phase approximation. We demonstrate that the 
correlation contribution to the compressibility in BLG is 
crucial when dielectric screening is weak and interactions 
within the graphene sheet are strong. Indeed, neglect of 
correlation effects leads to an error of the order of 100% 
in the case of suspended bilayers. We compare our re- 
sults for the compressibility of BLG with those obtained 
earlier for SLG'^" and are able to clearly identify the phys- 
ical origin of the main differences. For simplicity we as- 
sume here that the bilayer remains in a normal Fermi- 
liquid state down to very low densities. The behavior 
of the compressibility when one of the exotic low-density 
phases predicted in Refs. 21-25 is approached from the 
high-density Fermi-liquid phase is beyond the scope of 
the present theory. 

We note that the compressibility of BLG has already 
been calculated at the Hartree-Fock (HF) level in Ref. 13. 
We comment on the relationship between our results and 
those obtained in this earlier work in Sect. IIIB. We 
restrict our attention in this article to the case of a bal- 
anced bilayer in which inversion symmetry is not broken 
by an electrical potential difference between the layers. 
When a potential difference is present a gap opens up 
in the single-particle energy spectrum^; there is no gap 
between conduction and valence bands in the balanced 
bilayer limit that we consider. We also neglect trigonal 
warping effects in the bands which become important 
only at very low densities at which disorder effects nor- 
mally dominate. Both limitations are shared with the 
HF theory of Ref. 13. 

Our paper is organized as follows. In Sect. II we 
introduce the model and the linear-response functions 



which control BLG ground-state properties. In Sect. Ill 
we (i) derive expHcit expressions for the exchange and 
correlation energies using the integration-ovcr-coupling- 
constant algorithm and the fluctuation-dissipation theo- 
rem, (ii) introduce the random phase approximation for 
the correlation energy, and (iii) present and comment on 
our main numerical results. Finally, in Sect. IV we sum- 
marize our main findings and discuss their signifigance. 
Some technical details are relegated to an appendix. 



II. MODEL HAMILTONIAN AND 
LINEAR-RESPONSE FUNCTIONS 

BLG is modeled as two SLG systems separated by a 
distance d and coupled by both inter-layer hopping and 
Coulomb interactions. Most of the properties we dis- 
cuss below depend qualitatively on the Bernal stacking 
arrangement in which one sublattice (say A) of the top 
layer is a near-neighbor of the opposite sublattice (say 
B) of the bottom layer. Neglecting trigonal warping, the 
continuum model single-particle Hamiltonian^'^^ of a sin- 
gle valley is {h= 1), 



k,a,l3 



where Tap{k) are the coefficients of the following 4x4 
matrix 

r(fc) = -w7V7-fe-— (7''^7''+*7^) • (2) 



Here v (~ 10^ m/s) is the Fermi velocity of an isolated 
graphene layer, t± (~ 0.35 eV) is the inter-layer hopping 
amplitude, and the 7'' are 4x4 Dirac 7 matrices in the 
chiral representation'^^ (7^ = — i7°7^7^7^). The Greek 
indices a, j3 account for the sublattice degrees of freedom 
in top (1 = A, 2 = B) and bottom (3 = A, 4 = _B) layers. 
In the other valley the kinetic Hamiltonian is given by 

r'(fe) = r*(-fc).' 

If BLG is embedded in a medium with uniform dielec- 
tric constant e, electrons in the same layer interact via the 
2D Coulomb potential Va,{q) = 27re^/e(7, while electrons 
in different layers interact via VD(g) = Vs((7) exp (— qd). 
If the dielectric media above, below, and between the 
two layers are not identical these simple expressions are 
no longer valid"^^. 

For practical calculations of thermodynamic quanti- 
ties and linear-response functions it is convenient to 
work in the single-particle Hamiltonian eigenstate basis. 
Diagonalization of T(fc) yields four hyperbolic bands^^ 
with dispersions, ei,2(fc) = ±yjv'^k'^ + 1"]_/4 + t^_/2 and 
e3 4(fc) = ±-\/w2fc2 -\-t]_/A — t^/2. The interaction con- 
tribution to the Hamiltonian is 



T~L\nt — 



25^ [ 



V+{q)Pqp-q + V7(q)TqT_q 



(3) 



where S is the 2D electron system area, V± — (Vs ± 
Vd)/2, and pq and Tq are respectively the operators for 
the sum and difference of the individual layer densities: 



Pq 



= E 4-,,A("L,"fe) 



AA'Cfc.A' 



fe,A,A' 



and 



T — 



E 

k,X,X' 



k — q.X 



{1^l-ql^l^k)\yCk,y 



(4) 



(5) 



Here lAk is the unitary transformation from sublattice to 
band labels A, A' (see Appendix A). 

We evaluate interaction energies using a coupling- 
constant-integration scheme which expresses energies 
in terms of electronic equal-time correlation functions. 
The correlation functions can then be related to re- 
sponse functions using the fluctuation-dissipation theo- 
rem. When this commonly used approach^^ is adapted 
to the case of BLG, we see from "Hint that two response 
functions are necessary for the evaluation of ground-state 
properties of BLG^'^: the total-density response function. 



X+(.<1-,^) == -g{{Pq'-,P-q))oJ , 



\^) and the density-difference response function 



X-{q,^) = ^((Tq;T_q))„ . 



(6) 



(7) 
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Here {{A; 13))^ is the Kubo product 

p+oo 

{{A;B))^ = -t \im / dte'^'e-^'{[A{t),B]) , (8) 

(...) being the ground-state expectation value. 

At this point the reader might wonder why we have 
not introduced the mixed sum and difference response 
functions, ((/5q; T_q))„/5 and ((Tq;p_q))^/5'. As ex- 
plained in Ref. 19, these response functions vanish be- 
cause the system Hamiltonian is invariant under spatial 
inversion (parity). This is easily seen in the sublattice 
and layer basis where the parity operator V is given 
by P = (7^7'')*: with * indicating complex conjuga- 
tion. Using this compact expression for the parity op- 
erator "P, we can conveniently calculate its effect on one- 



body operators, like fiq 

for example, in the following manner: 



Efc,a,/3 4-q,a-^"/3(fc, g)cfc,/3, 



Z-^k.a 



,13 ^k-q,a 



VaqV 



'"f''TA{k, q)j''j^]*^I^Ckj- In this way, it is 



easy confirm that the density-sum operator is even under 
parity, VpqV = pq, while the density-difference operator 
is odd, VTqV = -Tq. 

Consider now a mixed response function such as 
{{pq]T^q))i^. We can write it in the exact-eigenstate 
(Lehmann) representation^^ as 



((pq;t-q)).. = E 



p - p 



IT] 



(*m|/5q|*„) 



X (*„|T_q|*, 



(9) 



Here w„ 



En — Em are excitation energies, P„ 



exjp(—/3En)/Z [with /3 = (kBT) and Z the canoni- 
cal partition function] are Boltzmann factors, Onm = 
(^„|0|^m) arc matrix elements of the operator O, and 
the limit 77 ^ 0+ is understood. We now use that the ex- 
act eigenstates |^„) of the system Hamiltonian are also 
eigenstates of the parity operator since the the Hamilto- 
nian is parity invariant: V\'^n) — ±|^n)- We find that 



((p.;t-,». = J2 



P — P 



7n,n 

t-^ P — P 

\ ^ ^ rn -* 7 



irj 






irj 



{'^m\VpqV\^n) 



(*m|Pq|*^ 



X (*„|T_q|*„) 
= -((Pq;t_q», 



(10) 



where we have used that pq is even under parity while 
Tq is odd. It follows that {{pq-^T-q))^ = 0. 

In the next Section wc will use the two response func- 
tions x_|_(<7,cj) and x_(<7,a;) to calculate exchange and 
correlation contributions to the BLG equation of state, 
and thus to the compressibility. 



III. EXCHANGE AND CORRELATION 
CONTRIBUTIONS TO THE COMPRESSIBILITY 

A. Formal electron-gas theory 



The compressibility k is defined by 



27 



1 



,dp, n? d'^E 
dn S din? 



(11) 



where p, = dE/dN is the chemical potential of the inter- 
acting system, E is the total ground-state energy, and n 
is the total (electron) density"^*. 

Using the Hellman-Feynman coupling-constant- 
integration theorem^^ and the specific form of "Hint given 
above in Eq. (3) we find that the interaction contribution 
to the ground-state energy of BLG is given by 



-E'int — 



N 



E 



dX 



(2^ 



Vt{q) 



Sf\q)-l 



(12) 



;.(A) 



where S'j_ (q) are the even and odd parity electron static 
structure factors at coupling constant A. Appealing to 
the fluctuation-dissipation theorem we find that 



jW 



S'^'il) 



1 

7rn 



+00 



dn x''l\q,in) . (13) 



This form of the fluctuation-dissipation theorem takes 
advantage of the smooth behavior of the linear-response 
functions Xj_ ((7,ii7) along the imaginary axis, which 



simplifies the task of performing accurate numerical 
wavevector and frequency integrals. (Along this axis one 
does not have to worry about the the collective plasmon 
poles and subtle particle-hole continuum band-edge fea- 
tures which are present along the real-frequency axis.^^) 
Substituting Eq. (13) in Eq. (12) we obtain an expres- 
sion for the total ground-state energy of the interacting 
system: 



E = E. 



N 



1 

TTn 



E 

i=± 



+00 







dX 



d^q 



(27r)2 
dil xi^\q,tn) - I 



Vdq) 



(14) 



Eq being the trivial noninteracting kinetic energy. Fol- 
lowing the conventional procedures of electron-gas the- 
pj.y27,30^ we separate out the contribution that is first 
order in e^ (i.e. the "exchange" energy) by writing 
E — Eq + E^ + Ec- The exchange energy per electron, 
Ex "^ E^/N, is given by 



IE 



d'^q 



-I 



(277)2 



Viiq) 



-00 

dnx^^\,q,i^)- 



(15) 



where x± (<?; i^) are the response functions of the nonin- 
teracting system, which have been calculated in Rcf. 19. 
The correlation energy (per electron) , which by definition 
is the sum of all the terms of higher order in e^, is given 

by 



1 

277ri 



E 



dX 



(2^) 



+00 



dQ Tiiq^ i^) 



(16) 
where Ti{q,iQ) = Vi{q)Axi {q,ifl) and Axi {q,iQ) = 

Xi {q,in) — Xi (qii^)- Neglecting correlations {i.e. 
treating interactions to first order in e^) one obtains the 
HF result for the ground-state energy. ^^ 

In this article we treat correlations within the ran- 
dom phase approximation (RPA)^'' in which the response 
functions of the interacting system at coupling constant 
A are given by 



Xi 



''\q,^) = 



xf\q,^) 

l-XV,{q)xf\q.io) 



(17) 



After inserting Eq. (17) in Eq. (16) one can observe that 
the integration over the coupling constant A can be per- 
formed analytically with the result 



.RPA 



77-rj ^ ^ 



2im 



=±1 



d^q 



+00 



dn 



{v,{q)xf\q,^^) 



In 1 



[l~V,{q)x^;\q,tn)]} 



(18) 



Since the hyperbolic bands of BLG asymptotically be- 
come linear in momentum (and thus identical to those of 
SLG), the integrals over frequency fi in Eqs. (15) and (16) 
diverge. We thus proceed as in the single-layer case'^'^ and 
regularize the frequency integrals by subtracting from e^ 
and £^ the infinite (and physically irrelevant) contin- 
uum modle exchange and RPA correlation energies of the 
undoped system. In this way we introduce a regularized 
exchange energy: 



Se^ = - 



1 

2Trn 



E 

i=±i 



d^q 

(27r)2 



Ve{q) 



+ CXD 



dndxfH,q,m, 

(19) 
where we have introduced the regularized response func- 

tions, Sxf\q,^n) ^ xf\q,^n)-xf''\q,^^), xf''\q,^n) 
being the noninteracting response functions of the un- 
doped system^^. The regularized correlation energy is 
given by 



-RPA 



-E 



27rn 



d^q 



~\-oo 



dn 



In 



1-Veiq)xf\q,m 
l-Ve{q)x^^''Hq,tn)\ 



Vdq)Sxf\q,^^) 



(20) 



This expression can be used to evaluate changes in in- 
teraction energy with carrier density at densities that 
are small compared to inverse unit cell area values, and 
is therefore reliable in the density range of relevance to 
gated and doped BLG electronic systems. Eqs. (19) 
and (20) are the most important results of this work. 
Together with the results for the doped Xe (^j*^) f^^^d 
undoped Xi " (?: *^) dynamical response functions pre- 
sented in Ref. 19 they allow us to accurately evaluate the 
ground-state energy (per electron) of BLG and thus the 
compressibility. 

After the regularization procedure described above, 
the integrals over J7 in Eqs. (19) and (20) are finite 
but the ones over q diverge. These divergences must be 
regularized by introducing an ultraviolet cut-off gmax = 
£uiax/v- Below we choose Ef/v as the unit of momentum, 
where e-p is the Fermi energy: 



£f 



t_L 

2 



n, if £i(fc) is occupied 



-— -I- v'^nn, if ei{k) is empty 



(21) 



Thus the integrals over dimensionless wave vectors must 
be calculated up to a maximum value A = emax/ep, cor- 
responding to the highest energies at which the contin- 
uum model applies. We set £max to 
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FIG. 1: (Color online) Top panel: interaction energy per elec- 
tron (in eV) as a function of doping (in units of lO^'^ cm~'^) for 
different values of graphene's fine-structure constant Oee- The 
values of Ofoo displayed are (from bottom to top) Oee = 0.125 
(solid line), 0.25 (short-dashed line), 0.5 (dotted line), 1 
(long-dashed hue), 2.2 (dash-dotted line). Note the cusp for 
n~18xl0^^ cm~^, which is more prominent for large values 
of Cfeo, the value of doping at which the high-energy split off 
band ei(fc) is first occupied. Bottom panel: the chemical po- 
tential fi (in eV) as a function of doping for different values 
of Qcc. Color-coding and labeling is the same as in top panel. 
Crosses label n for Qcc ~ (noninteracting bilayer graphene) . 



£ m a.TC — 



/27rw2 



7.2 eV 



(22) 



3\/3aQ/2 is the SLG unit-cell area, oq 



where ^o 

1.42 A being the carbon-carbon distance, and n 



'i/Ao is the number of 7r-electrons per unit area in the 
neutral system. The energies we evaluate have a weak 
logarithmic dependence on £inax that has little impor- 
tance for the conclusions we draw below. 
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FIG. 2: (Color online) The Hartree-Fock (exchange-only) 
inverse thermodynamic density-of-states dfj,/dn (in units of 
eV X A ) as a function of the Fermi energy ep (in units of 
£max X 10"'') for different values of Qoc- The Fermi energy 
£f ranges from ep = 7 x 10~* eV (corresponding to a doping 
n ~ 2 X 10^" cm~'^) to ep = 0.53 eV (corresponding to a dop- 
ing n ~ 4.0 X 10^"' cm~^). Color-coding and labeling is the 
same as in Fig. 1. Crosses label d^/dn for Oee = (noninter- 
acting bilayer graphene). A negative (5-function contribution 
to dn/dn at n = Til for Oee 7^ has been omitted. 



B. Numerical results and discussion 

We now turn to our main numerical results. The 
ground-state properties of BLG are completely deter- 
mined by the total density n, by the interlayer distance 
d, which we have taken to be d = 3.35 A, by the inter- 
layer hopping t_L, which we have taken to be 0.35 eV, and 
by the fine-structure constant (restoring h for a moment) 

In the top panel of Fig. 1 we present the exchange- 
correlation energy 5e^c = <5ex + i5e|^^^ as a function of 
n and a^e- For our choice of energy zero, (5£x is positive 
and Sec is negative. The two contributions to the in- 
teraction energy tend to cancel strongly, with a slightly 
positive total, suggesting that correlation will be as im- 
portant as exchange in determining physical properties. 
We see that 5e^c has a cusp at every value of aoe for 
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FIG. 3: (Color online) Top panel: same as in Fig. 2 but 
with the inclusion of RPA correlations. The thick solid 
(cyan) line labels the inverse thermodynamic density-of-states 
of suspended (aeo = 2.2) single-layer graphene. The val- 
ues of ep range from ep = 7 x 10"** eV to ep = 0.53 eV 
(n « 4.0 X 10^^ cm~^). Note that the compressibility of bi- 
layer graphene remains finite for n — >■ 0, while the one of 
single-layer graphene diverges. A negative (5-function contri- 
bution to dfi/dn a,t n — ni for Ofoo 7^ has been omitted. 
Bottom panel: a zoom of the upper panel for low densities. 
Note that the horizontal axis in the bottom panel represents 



ni 



2(t_L/w)2/7r w 18 X 10^2 cm"2, the value of dop- total density n in units of 10 cm 



ing at which ep = t±- It is at this value of n that the high- 
energy split off band Si{k) is first occupied. As a conse- 
quence, the chemical potential ^ = (3[n((5£kin + <5exc)]/9n 
has a jump at n = ni when doping is increased from val- 
ues smaller than ni to values larger than ni. Here fekin 
is the kinetic energy (per electron) of the noninteracting 
system with the Dirac point chosen as energy zero. The 



chemical potential /i as a function of doping is illustrated 
in the bottom panel of Fig. 1. Note that the jump is 
downward. These type of jumps have been discussed ear- 
lier in the context of second-subband occupation in wide 



quantum wells (see for example Ref. 35 and references 
therein) and are potentially technologically interesting 
since they can in principle lead to bistability. We observe 
that the chemical potential jump implies a (5-function sin- 
gularity in the compressibility k at n = rii: this feature 
will be omitted in the presentation of compressibility nu- 
merical results below, Figs. (2)- (4). 

In Fig. 2 we report HF theory results for the inverse 
thermodynamic density-of-states dfi/dn which has ki- 
netic and exchange energy contributions: 9y^/9n|HF = 
9^[n((5£kin + '^£x)]/5n^. The decrease in d^ijdn with den- 
sity at aco = is a consequence of the difference be- 
tween hyperbolic and parabolic dispersion. We see that 
dji/dn is positive and enhanced by exchange interactions 
over the density range covered in this plot. The non- 
monotonic behavior predicted in Ref. 13 appears only at 
extremely low densities and, as we discuss below, does 
not survive RPA correlations. This behavior contrasts 
with that of ordinary 2D electron gases in which HF 
theory predicts a negative compressibility below a criti- 
cal density that is moderate, and only weakly influenced 
by correlations. This qualitative behavior difference is 
a consequence of the relevance of exchange interactions 
with both conduction and valence bands, and of the 
sublattice pseudospin chirality of the bands. The same 
mechanisms are also responsible for a thermodynamic 
density-of-states that is suppressed rather than enhanced 
in SLG^^'^°. Like an ordinary 2D electron gas, BLG has 
an intra-conduction-band exchange contribution to its 
chemical potential that is negative and proportional to 
n}'"^; this energy is however approximately half as large in 
the BLG case because the wavefunctions are spread over 
more than one sublattice. The tendency toward a nega- 
tive compressibility is further countered in the BLG case 
by inter-band exchange, which yields a chemical poten- 
tial contribution that is also proportional to v}/"^ in the 
low-density limit, but positive. The end result is that the 
low-carrier density negative ■n}''^ exchange contribution 
to the chemical potential, which would yield a negative 
compressibility, is approximately six times smaller (for 
the same background dielectric constant) in BLG than 
in an ordinary 2D electron gas. When only exchange in- 
teractions are included, we find that the total compress- 
ibility calculated within the two-band model^ becomes 
negative only at densities below a critical value given by 
(restoring h for a moment) 



1 ( I3accti_ 
2hv 



al^{l.l X 10^") cm" 



(23) 



in agreement with the numerical results in Fig. 2 of 
Ref. 13. In Eq. (23) we have introduced 



P- 



1 



3 
16 



1 2 

dx^ 2^1(5/2, 1/2, 3, l/x2^ 



9tt ' 
(24) 
2-^1 (a, 6, c, x) being the usual hypergcomctric function. 
For the sake of comparison, note that within HF the crit- 
ical density at which the compressibility of a standard 2D 



electron gas changes sign is nj; = 2/{TT^a%), where 

Ob is the material Bohr radius. In GaAs, for example, 
an « 100 A and thus ni^^^^^ w 6.5 x 10^" cm^^. In the 
graphene case we find numerically that in the random 
phase approximation the low-density negative compress- 
ibility does not survive correlations. This compressibil- 
ity anomaly is in any event likely to be preempted by 
BLG's low-density ferroelectric instability^^'^^'^"*, which 
is driven by physics beyond that captured by the RPA as 
discussed above. The issue of a possible negative com- 
pressibility in BLG is discussed further below. 

In Fig. 3 we report on results for the inverse ther- 
modynamic density-of-states, dix/dn, calculated includ- 
ing both exchange and RPA correlations corrections: 
dfi/dn = 9^[n(fckin + Se^c)]/dn'^. Qualitatively, the re- 
sults in Fig. 3 look rather similar to the HF ones in Fig. 2. 

However, we clearly see that for dopings below ni « 
18 X 10^^ cm^^ correlation effects are quantitatively very 
important. For instance, percentage values of the ratio 



r{acc) = lim 



d^/dn 



{dfi/dn)\ 



HF 



(25) 



(between the data in Fig. 3 and the HF data in Fig. 2) 
are of the order of « 20% for a^c — 0.5 and larger than 
100% for a suspended bilaycr (ctoc = 2.2). 

For the sake of comparison in Fig. 3 we have also plot- 
ted dfi/dn for suspended (ofoc = 2.2) SLG. As expected, 
the difference between DLG and SLG compressibilities 
is very small at high doping, especially so when all four 
bands Si{k) are occupied. At low densities, however, the 
results are very different since in this regime the BLG 
spectrum approaches a parabolic form k'^/{2m), with 
m = t±/(2v'^), strongly deviating from the SLG linear 
dispersion. In particular, note that d^/dn diverges for 
n — > in the SLG case, while it approaches a finite value 
for BLG. This striking difference stems from the behav- 
ior of the BLG quasiparticle effective mass to*, which 
remains finite when doping approaches zero^^. 

In Fig. 4 we plot the ratio k/kqj '^o = 
[n?d'^{n5e\^i-a)/dii?']^^ being the compressibility of the 
noninteracting system. We clearly see from this plot that 
the main effect of electron-electron interactions is to sup- 
press n. This can be easily understood within the Landau 
theory of normal Fermi liquids. Indeed k/kq is largely 
controlled by and proportional to m* /m and, as demon- 
strated in Ref. 16, the role of interactions is to suppress 
TO.* with respect to the bare value, i.e. m* jm < 1. As ex- 
plained in Ref. 30 for the case of SLG, the suppression of 
the mass (or enhancement of the quasiparticle velocity) 
stems from the chiral nature of the low-energy spectrum. 

We now turn to a discussion of the compressibility 
in the extreme low-density limit, illustrated in Fig. 5. 
As discussed above and first explained by Kusminskiy et 
al.^'^, the compressibility becomes negative at extremely 
low densities in the HF approximation because of a small 
net contribution to the chemical potential that is negative 
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FIG. 4: (Color online) The ratio between the interacting- 
system compressibility (k) and the noninteracting one (kq), 
as a lunction of doping (in units of lO"'^^ cm"^). Color coding 
and labeling are the same as in Figs. 2-3. Note that the ratio 
k/kq is smaller than unity, the more so the stronger electron- 
electron interactions are. 



and proportional to n^". This contribution to the com- 
pressibility is reminiscent of the larger but related con- 
tribution to the chemical potential which appears in an 
ordinary 2D electron gas. Because the relative strength 
of interactions and band-energies in that case can be ab- 
sorbed in a length scale change, the n^/^ exchange energy 
can be viewed as the leading order term in an expansion 
of energy in powers of e^/fcp ^ e^n~^". The leading or- 
der correlation contribution to the chemical potential is 
therefore proportional to n'^ (up to logarithmic factors) 
and does not appear in the compressibility. This sim- 
ple scaling property does not apply to BLG, because the 
inter-layer hopping and in-plane hopping terms in the 
continuum-limit Hamiltonian do not scale in the same 
way with density. The chemical potential and energy per 
particle have to be expanded separately in terms of pow- 
ers of n^/^ and the interaction scale a^e- As illustrated 
in Fig. 5 we find numerically that both the exchange 
and correlation contributions to the chemical potential 
change sign at very low carrier densities, in such a way 
that the total chemical potential is a monotonically in- 
creasing function of energy. In fact we find that the ratio 
k/kq is smallest at low density: in other words, as it is 
also clear from Figs. 2-3, the enhancement of dfi/dn rel- 
ative to the noninteracting system results becomes larger 
for lower carrier densities. 



The competing role of exchange and correlations here 
is similar to their role in the spin-susceptibility. The ex- 
change energy favors spin-polarization by lowering the 
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FIG. 5: (Color online) The various contributions to the 
ground-state energy (per particle and divided by n^") as 
functions of the carrier density n in units of 10^^ cm~^. These 
results refer to Oee = 1. 



chemical potential of each spin when its occupation in- 
creases. When correlations are included, the energy de- 
pends more on the total density and less on its partition- 
ing into spin components. Similarly here the exchange 
energy favors occupation of the higher subband at low- 
densities because of the large Coulomb interaction matrix 
elements when the Fermi radius is small. When corre- 
lations - which are less sensitive to interactions between 
quasiparticles at the Fermi energy - are restored, the sen- 
sitivity to band index is reduced and the overall trend in 
the dependence of chemical potential on density is re- 
stored. 



IV. SUMMARY AND DISCUSSION 

In summary, we have calculated the compressibil- 
ity of crystalline bilayer graphene beyond the Hartree- 
Fock approximation by treating correlation effects at the 
random-phasc-approximation level. We have shown that 
electron-electron interactions suppress the compressibil- 
ity quite substantially with correlation effects playing an 
important quantitative role. The reduction of compress- 
ibility stands in stark contrast to the large compressibil- 
ity enhancements that occur in regular two-dimensional 
electron gas systems, even though the two systems share 
the same parabolic dispersions. The source of the qual- 
itatively different behavior is the importance in bilayer 
graphene of exchange interactions between carriers in the 
conduction band and the full negative energy Dirac sea. 
The suppression of the compressibility has the same ori- 
gin as the enhancement of quasiparticle effective mass. 
Both phenomena ultimately originate from the chiral na- 



ture of the low-energy spectrum. 

The present results demonstrate that correlations play 
an essential role in quantitative studies of interaction 
effects in bilayer graphene. Previous work^^'^^'^^ has 
suggested that neutral bilayers might become unsta- 
ble to spontaneous layer polarization when disorder is 
weak. The role of long-range Coulomb interactions in the 
physics of this instability could be addressed by extend- 
ing the calculations described here to the case in which 
there is an electric potential difference between layers. 
One complication associated with this elaboration is that 
inversion symmetry is explicitly broken so that correla- 
tions between even and odd parity density fluctuations 
are non-zero. The more complicated form for the sub- 
band spinors of the band eigenstates also makes the task 
of finding quasi-analytic results for the noninteracting 
system polarization (Lindhard) functions challenging. In 
the present calculation the semi- analytic results we have 
used for the Lindhard function^^ are extremely helpful 
and allow the wavevector and frequency integrals to be 
evaluated numerically with confidence and precision, in 
spite of the numerical subtleties that lurk in the inte- 
grands. It is difficult to obtain accurate results when the 
Lindhard function is evaluated by brute-force numerics. 
Although this lies outside the scope of the present paper, 
we anticipate that correlation effects will lower predic- 
tions for the amount of charge transferred between the 
layers. 

After this work was complete we learned of three recent 
experimental studies'^^^'^^ which measure the compress- 
ibility of BLG. All three groups find that, in the bal- 
anced limit, dfi/dn has a peak near zero carrier density 
and then decreases monotonically with increasing car- 
rier density; the change of sign at low densities which 
appears in an ordinary two-dimensional electron gas is 
absent in BLG, in agreement with our findings. Our 
calculations demonstrate that these experimental results 
can be strongly influenced by interactions so that some 
caution must be exercised in fitting compressibility mea- 
surements to band-structure models. 
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Appendix A: The unitary transformation Uh 

For the sake of completeness, in this Appendix we re- 
port the form of the unitary matrix Uk which diagonalizes 
the kinetic matrix T{k). The matrix Uk can be written 
as 



Uk = G{cf,k)AR{0k) 



(Al) 



where the two angles (j)k and Op, are defined by (j)^ ~ 
a,TctSin{ky / kx) and 9k = arctan[Ay£3(fc)/ei(fc)]. 

The 4x4 matrices G{(t)k), A, and R{9k) in Eq. (Al) 
are given by 



Gal3{4>k) = 5al3[5al + S^A - 6 '5a2 - e *'*''°(5a3] , (A2) 



-i<t>ki 



dij being the usual Kronecker delta. 



1 



^=^(7'+7V), 



(A3) 



and, finally, R{9k) is the tensor product of two 2 x 
2 rotations which act on the subspace spanned by 
the odd (inversion-antisymmetric) and even (inversion- 
symmetric) eigenstates of the kinetic Hamiltonian T, re- 
spectively. Recalling that the odd bands arc those corre- 
sponding to the eigenvalues £1,2(^)7 while the even ones 
are those labeled by the eigenvalues £3^4 (fc), we find that 
R{9k) = R2{9k)\i,2 ® R2{9k%A with ' 
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cos{9k) sin(6'fc) 
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